Baf53cre ENS neurons, SI data
Steady state, CTL x3 CKO(Ahr-cko) x3
loading 60k cells,
cellranger called 24,365 cells (demo); 24,604 cells (plus)
pure neurons downstream
GEX.seur <- readRDS("./GEX0429.seur.preAnno.rds")
GEX.seur
## An object of class Seurat
## 22951 features across 13503 samples within 1 assay
## Active assay: RNA (22951 features, 1040 variable features)
## 3 dimensional reductions calculated: pca, tsne, umap
color.FB <- ggsci::pal_d3("category20c")(20)[c(1,6,11,#16,
2,7,12#,17
)]
color.cnt <- color.FB[c(1,4)]
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno1") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
GEX.seur <- subset(GEX.seur, subset = preAnno1 %in% levels(GEX.seur$preAnno1)[1:16] & DoubletFinder0.1 == "Singlet")
GEX.seur
## An object of class Seurat
## 22951 features across 11932 samples within 1 assay
## Active assay: RNA (22951 features, 1040 variable features)
## 3 dimensional reductions calculated: pca, tsne, umap
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno1") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
markers.new.s <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
"Gfra2","Oprk1","Adamtsl1",
"Fbxw15","Fbxw24","Chrna7",
"Satb1","Cntnap5b",
"Gabrb1","Nxph1","Lama2","Lrrc7",
"Ryr3","Eda","Tac1",
"Kctd8","Ntrk2","Penk",
"Fut9","Nfatc1","Egfr","Mgll",
"Chrm3"
),
IMN=c("Nos1","Kcnab1",
"Gfra1","Etv1",
"Man1a","Airn",
"Adcy2","Cmah","Creb5","Vip","Pde1a",
"Ebf1","Gpc5"
),
IN=c("Npas3","Synpr","St18","Gal",
"Kcnk13",
"Neurod6","Moxd1","Sctr",
"Piezo1","Sst","Adamts9","Kcnn2",
"Calb2"),
IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
"Ngfr","Galr1","Il7","Aff2",
"Gpr149","Cdh6","Cdh8",
"Clstn2","Ano2","Ntrk3",
"Cpne4","Vwc2l","Cdh9","Scgn",
"Vcan","Cck","Piezo2","Kcnh7",
"Rerg","Bmpr1b","Skap1","Ntng1",
"Tafa2","Nxph2"),
CONTAM=c("Actl6b","Phox2b"),
pDEGs=c("Grid2","Ide","Btaf1","Cdh10","Slc15a2","Ccser1","Hdac9","Col25a1","Rspo2"))
pm.All_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new.s
pm.CTL_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CTL")), features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CTL only")
pm.CTL_new.s
pm.CKO_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("CKO")), features = as.vector(unlist(markers.new.s)), group.by = "preAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CKO only")
pm.CKO_new.s
GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1500)
# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2
head(VariableFeatures(GEX.seur), 200)
## [1] "Mgat4c" "Cntn5" "Zfp804a" "Klhl1"
## [5] "Wdr17" "Gm30613" "Cntnap2" "Kcnb2"
## [9] "Adgrg6" "Robo2" "Gal" "Brinp3"
## [13] "Cpne4" "Lingo2" "Cdh8" "Zfp804b"
## [17] "Nmu" "Nrxn3" "Ntng1" "Kcnip4"
## [21] "Gm21847" "Gm32647" "Gm29521" "Ano2"
## [25] "Pcdh9" "Gm15261" "Pdzrn4" "Grm5"
## [29] "Sgcz" "Tmeff2" "Prkg1" "Ntrk3"
## [33] "Cdh9" "Gpr149" "AC150683.1" "Nrg1"
## [37] "Ebf1" "Cdh18" "Schip1" "Dgkg"
## [41] "Car10" "Cmah" "Kcnq5" "Fgf14"
## [45] "March1" "Astn2" "4930509J09Rik" "Vwc2l"
## [49] "Cdh6" "Nxph2" "Clstn2" "Vip"
## [53] "Asic2" "Sst" "Dcc" "Sema5a"
## [57] "Gm20754" "Unc5d" "Grin3a" "Spock3"
## [61] "Gm26612" "Piezo2" "4930432L08Rik" "Dpp10"
## [65] "Gpc5" "Gm15680" "1700063D05Rik" "Efr3a"
## [69] "Pbx3" "Pcdh10" "Egfem1" "Serpini1"
## [73] "Csmd3" "Gm26737" "Fut9" "Ccbe1"
## [77] "Kcnh7" "Zfhx3" "Gm15584" "Itgb6"
## [81] "Gm43388" "Rasgef1b" "Gm32916" "Dgki"
## [85] "Dgkb" "4930422I22Rik" "Plxna4" "Myo3b"
## [89] "Cysltr2" "Plcl1" "Ttc29" "Tbx20"
## [93] "1700034P13Rik" "Chrna9" "1700111E14Rik" "9530059O14Rik"
## [97] "Tnr" "Chsy3" "Cpa6" "Pde1a"
## [101] "Pcdh11x" "Robo1" "Gm20713" "Gm49938"
## [105] "D030068K23Rik" "Penk" "Luzp2" "B230323A14Rik"
## [109] "Pde4d" "Csmd1" "Il1rapl1" "Ssc4d"
## [113] "Col25a1" "Trps1" "Casp8" "Trhde"
## [117] "Myl1" "Abca9" "Arhgap6" "Cacna2d3"
## [121] "Grik1" "Ghr" "Rab27b" "5530401A14Rik"
## [125] "Galntl6" "Tafa2" "Calcb" "Gm11342"
## [129] "B230312C02Rik" "Skap1" "Gm19784" "Kcnd2"
## [133] "4931415C17Rik" "Tac1" "Nwd1" "Gm12239"
## [137] "Galr1" "Col8a1" "Ccdc3" "Gm28494"
## [141] "Clcnka" "Hsd17b14" "Hs6st3" "Gna14"
## [145] "4930471E19Rik" "Cntn4" "Chrm3" "Adamts9"
## [149] "Cux2" "Kctd8" "Lama2" "Gm20560"
## [153] "Gm26713" "Gm50255" "Gm37267" "Galnt18"
## [157] "Slc4a4" "Nos1" "Nxph1" "4930587E11Rik"
## [161] "4930445B16Rik" "A330008L17Rik" "Gm49127" "Tenm2"
## [165] "Rarb" "Gm49678" "Aldh1a2" "Cdh10"
## [169] "2310039L15Rik" "Cadm2" "Kif26b" "Gabrg3"
## [173] "Gm30094" "Zfp286os" "Gm11339" "Pdgfra"
## [177] "Gm16685" "Nox3" "Bmpr1b" "Apol7b"
## [181] "Dock10" "Adamts6" "Gm34544" "Fam81b"
## [185] "Adam2" "Gm49953" "Cbln2" "Mid1"
## [189] "Ccnb1" "Npas3" "Cntn3" "C79798"
## [193] "Kctd16" "Atp12a" "Gm16226" "Tafa1"
## [197] "Synpr" "Cep72" "Adap2os" "Gm38405"
GEX.seur <- ScaleData(GEX.seur, features = rownames(GEX.seur))
## Centering and scaling data matrix
# exclude MT genes and more
#DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^Gm|^Hist",rownames(GEX.seur),value = T)
DIG <- grep("^mt-|^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^AC|^AA|^AI|^AW|^BC|^Gm|^Hist|Rik$|-ps",
rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
c(#MT_gene,
DIG,
CC_gene) )
GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur))
## PC_ 1
## Positive: Lrrtm4, Galntl6, Tox, Ndst4, Epha5, Kcnd2, Grik1, Bnc2, Zfp521, Pcdh7
## Nos1, Pde3a, Kcnab1, Lrp1b, Cdc14a, Tshz2, Chrna7, Zfp536, Syt6, Specc1
## Synpo2, Plcxd3, Dach1, Lrrc4c, Fbxw15, Rspo2, Brinp2, Fat3, Tenm3, Fbxw24
## Negative: Ano2, Ntrk3, Robo2, Tmeff2, Cdh8, Cpne4, Nrxn3, Adgrg6, Gpr149, Pdzrn4
## Clstn2, Cntn5, Mgat4c, Spock3, Cdh6, Zfp804a, Dgkg, Grin3a, Pcdh10, Sgcz
## Cysltr2, Myl1, Arhgap6, Astn2, Plxna4, Cux2, Cacna2d3, Itgb6, Iqgap2, Tnr
## PC_ 2
## Positive: Rbfox1, Bnc2, Ptprt, Tshz2, Tafa1, Gpc6, Grik1, Frmd4b, St6galnac3, Brinp2
## Fbxw15, Caln1, Fbxw24, Pcdh7, Tox, Plcxd3, Tmem132c, Cdc14a, Tcf7l2, Specc1
## Pld5, Unc5c, Dock2, Adamtsl1, Oprk1, Alk, Sdk2, Nfia, Tmem132cos, Ccbe1
## Negative: Nos1, Auts2, Egfem1, Schip1, Kcnq5, Asic2, Cadps2, Etv1, Kcnab1, Epha5
## Cmah, Alcam, Dgkb, Gfra1, Dach1, Kcnt2, Hs6st3, Rgs6, Ebf1, Lrrc4c
## Stxbp6, Cntnap5a, Il1rapl1, Rgs7, Kcnj3, Cdh11, Pde4d, Fat3, Creb5, Kitl
## PC_ 3
## Positive: Cdh18, Klhl1, Kcnip4, Pbx3, Kctd8, Csmd3, Gpc6, Htr4, Cadm2, C79798
## Meis1, March1, Vwc2l, Sema5a, Zfhx3, Dlc1, Gabrg3, Serpini1, Zbbx, Galnt18
## Cdh9, Cntn3, Skap1, Csmd1, Car10, Plcl1, Khdrbs2, Tenm2, Stxbp5l, Sybu
## Negative: Adgrg6, Sgcd, Nos1, Nmu, Cysltr2, Gpr149, Rgs6, Grin3a, Slc4a4, Cntnap5a
## Itgb6, Kcnab1, Nfia, Dapk2, Ano2, Fgf13, Cpne4, Ngfr, Ccbe1, Cbln2
## Zfp804a, Efr3a, Dgkg, Zfp536, Gfra1, Otof, Syt2, Lhfpl2, Calcb, Pcdh10
## PC_ 4
## Positive: March1, Klhl1, Sema5a, Chsy3, Vwc2l, Cdh9, Dcc, Kcnh7, Galnt18, Ebf1
## Ntng1, Prune2, Kcnb2, Zfhx3, Pbx3, Bcl2, Rasgef1b, Tenm4, Trhde, Cntn4
## Zbbx, Il1rapl1, Sdk1, C79798, Alk, Cpa6, Gal, Sez6l, Col18a1, Lncbate1
## Negative: Dock10, Lingo2, Kcnip4, Ndst4, Tac1, Gda, Prkg1, Sorcs1, Nxph1, Epha5
## Dmd, Thsd4, Kcnt2, Csmd3, Lrrtm4, Cntn5, Lrrc7, Kctd8, Lrrc4c, Lrp1b
## Chl1, Cntn3, Fgf13, Lin7a, Lama2, Pgm5, Hgf, Galntl6, Pld5, Ccdc60
## PC_ 5
## Positive: Trhde, Cntn4, Ptprd, Trps1, Nrg1, Ebf1, Cpa6, Csmd1, Kcnd2, Gal
## Col18a1, Zmat4, Npas3, Egfem1, Rmst, Luzp2, Kcnk13, Fstl4, Csmd3, Nkain3
## Nav2, Chsy3, Synpr, Hs6st3, Cntn3, Shisa6, Adgrl2, Myo1e, Lrp1b, Hcn1
## Negative: Dgkb, Alcam, Kcnt2, Klhl1, Thsd7b, Vwc2l, Dpp10, Il1rapl1, Cdh9, Mgat4c
## Rasgef1b, Pbx3, Fgf13, Nos1, Lrrc4c, Zfhx3, Galnt18, Auts2, Khdrbs2, Slc44a5
## Zbbx, C79798, Vcan, Olfr78, Rgs6, Stxbp6, Nmur2, Cdh20, P3h2, Hmcn1
length(setdiff(VariableFeatures(object = GEX.seur),
c(#MT_gene,
DIG,CC_gene)))
## [1] 1027
setdiff(VariableFeatures(object = GEX.seur),
c(#MT_gene,
DIG,CC_gene))[1:300]
## [1] "Mgat4c" "Cntn5" "Zfp804a" "Klhl1" "Wdr17"
## [6] "Cntnap2" "Kcnb2" "Adgrg6" "Robo2" "Gal"
## [11] "Brinp3" "Cpne4" "Lingo2" "Cdh8" "Zfp804b"
## [16] "Nmu" "Nrxn3" "Ntng1" "Kcnip4" "Ano2"
## [21] "Pcdh9" "Pdzrn4" "Grm5" "Sgcz" "Tmeff2"
## [26] "Prkg1" "Ntrk3" "Cdh9" "Gpr149" "Nrg1"
## [31] "Ebf1" "Cdh18" "Schip1" "Dgkg" "Car10"
## [36] "Cmah" "Kcnq5" "Fgf14" "March1" "Astn2"
## [41] "Vwc2l" "Cdh6" "Nxph2" "Clstn2" "Vip"
## [46] "Asic2" "Sst" "Dcc" "Sema5a" "Unc5d"
## [51] "Grin3a" "Spock3" "Piezo2" "Dpp10" "Gpc5"
## [56] "Efr3a" "Pbx3" "Pcdh10" "Egfem1" "Serpini1"
## [61] "Csmd3" "Fut9" "Ccbe1" "Kcnh7" "Zfhx3"
## [66] "Itgb6" "Rasgef1b" "Dgki" "Dgkb" "Plxna4"
## [71] "Myo3b" "Cysltr2" "Plcl1" "Ttc29" "Tbx20"
## [76] "Chrna9" "Tnr" "Chsy3" "Cpa6" "Pde1a"
## [81] "Pcdh11x" "Robo1" "Penk" "Luzp2" "Pde4d"
## [86] "Csmd1" "Il1rapl1" "Ssc4d" "Col25a1" "Trps1"
## [91] "Casp8" "Trhde" "Myl1" "Abca9" "Arhgap6"
## [96] "Cacna2d3" "Grik1" "Ghr" "Rab27b" "Galntl6"
## [101] "Tafa2" "Calcb" "Skap1" "Kcnd2" "Tac1"
## [106] "Nwd1" "Galr1" "Col8a1" "Ccdc3" "Clcnka"
## [111] "Hsd17b14" "Hs6st3" "Gna14" "Cntn4" "Chrm3"
## [116] "Adamts9" "Cux2" "Kctd8" "Lama2" "Galnt18"
## [121] "Slc4a4" "Nos1" "Nxph1" "Tenm2" "Rarb"
## [126] "Aldh1a2" "Cdh10" "Cadm2" "Kif26b" "Gabrg3"
## [131] "Zfp286os" "Pdgfra" "Nox3" "Bmpr1b" "Apol7b"
## [136] "Dock10" "Adamts6" "Fam81b" "Adam2" "Cbln2"
## [141] "Mid1" "Ccnb1" "Npas3" "Cntn3" "C79798"
## [146] "Kctd16" "Atp12a" "Tafa1" "Synpr" "Cep72"
## [151] "Adap2os" "Loxl1" "Galnt13" "Met" "Ano5"
## [156] "Otof" "Ngfr" "Alcam" "Rbfox1" "Flrt2"
## [161] "Lrrc4c" "Grm7" "Zbbx" "Bmp5" "Hccs"
## [166] "Mme" "Scnn1a" "Hcn1" "Acp7" "Iqgap2"
## [171] "Cadps2" "Thsd4" "Plxdc2" "Eya1" "Fbxl7"
## [176] "Thsd7b" "Grp" "C1qtnf7" "Antxr2" "St8sia2"
## [181] "Tenm4" "Hsd17b3" "Foxa3" "Fmo2" "Slc13a1"
## [186] "Gpc6" "Ptprt" "Auts2" "Ggt5" "Bnc2"
## [191] "Nell1" "Calcrl" "Pigk" "Necab1" "Ptprd"
## [196] "Agtr1b" "Dapk2" "Olfr78" "Rasgrf1" "Prlr"
## [201] "Fam47e" "Kng1" "Sorcs1" "Hormad1" "Ptk7"
## [206] "Wee2" "Dach1" "Rerg" "Cdhr1" "Cpne8"
## [211] "Hs3st4" "Scn7a" "B3galt1" "Col18a1" "Bcl2"
## [216] "Tnfrsf23" "Tenm1" "Tex14" "Hnf1aos1" "Serpine2"
## [221] "Tex21" "Htr4" "Slc2a13" "Spag16" "Serpinb5"
## [226] "Qrfpr" "Gad2" "Cavin2" "Aff2" "Sorcs3"
## [231] "Naaladl2" "Prr16" "Agrn" "Esr1" "St6galnac5"
## [236] "Mgat4e" "Nmur2" "Nell1os" "Chst9" "Prune2"
## [241] "Slc22a20" "Lhfpl2" "Muc16" "Tgfb1" "Ror1"
## [246] "Stxbp6" "Rgs6" "Lrp1b" "Meis2" "Cemip"
## [251] "Zbtb7c" "Ppp3ca" "Pbld1" "Mlph" "Kif28"
## [256] "Wt1os" "Ovol2" "Cga" "Vmn2r18" "Akr1d1"
## [261] "Arhgap25" "Adprhl1" "Itgb2" "Pnoc" "Mast4"
## [266] "Rmst" "Prkca" "Etv1" "Vcan" "Caln1"
## [271] "Lmcd1" "Hs3st2" "Rtl4" "Tmem132d" "Plcb1"
## [276] "Yae1d1" "Sctr" "Brinp2" "Dmd" "Gfra1"
## [281] "Moxd1" "St18" "Adcy2" "Nkain3" "Adamtsl5"
## [286] "Slc25a41" "Fhit" "Pla2g10" "Epha8" "Dach2"
## [291] "Ndst4" "Ildr2" "Robo3" "Thpo" "Adgrl2"
## [296] "Fbxw24" "Sybu" "Shisa9" "Pde7b" "Cntnap5b"
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB) +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB)
DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)
ElbowPlot(GEX.seur,ndims = 50)
PCs <- 1:24
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 2)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11932
## Number of edges: 520184
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8052
## Number of communities: 26
## Elapsed time: 1 seconds
GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 145)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 00:18:27 UMAP embedding parameters a = 0.9922 b = 1.112
## 00:18:27 Read 11932 rows and found 24 numeric columns
## 00:18:27 Using Annoy for neighbor search, n_neighbors = 20
## 00:18:27 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 00:18:28 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\RtmpGCQCJl\file41804e884ea1
## 00:18:28 Searching Annoy index using 1 thread, search_k = 2000
## 00:18:31 Annoy recall = 100%
## 00:18:31 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 00:18:32 Initializing from normalized Laplacian + noise (using irlba)
## 00:18:33 Commencing optimization for 200 epochs, with 355450 positive edges
## 00:18:46 Optimization finished
DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno2") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))
length(as.vector(unlist(markers.new.s)))
## [1] 92
FeaturePlot(subset(GEX.seur, subset = cnt %in% c("CTL")),
features = as.vector(unlist(markers.new.s)), cols = c("lightgrey","red"))
FeaturePlot(subset(GEX.seur, subset = cnt %in% c("CKO")),
features = as.vector(unlist(markers.new.s)), cols = c("lightgrey","red"))
GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
levels = c(0,1,2,12,11,17,9,10,18,
4,5,14,8,7,13,16,23,
6,19,24,
3,20,
15,21,
25,22))
table(GEX.seur@meta.data[,c("preAnno2","sort_clusters")])
## sort_clusters
## preAnno2 0 1 2 12 11 17 9 10 18 4 5 14 8 7
## EMN1 1197 856 706 420 295 0 5 0 0 0 0 0 0 0
## EMN2 0 0 0 0 145 188 7 0 0 0 0 0 0 0
## EMN3 0 0 0 0 2 137 503 41 4 0 0 0 0 0
## EMN4 0 0 0 0 0 1 13 482 4 0 0 0 0 0
## EMN5 0 0 0 0 1 0 1 2 308 0 0 0 0 0
## IMN1 0 0 0 0 0 0 0 0 0 583 564 311 542 455
## IMN2 0 0 0 0 0 0 0 0 0 0 0 97 1 90
## IMN3 0 0 0 0 0 0 0 0 0 0 1 0 0 0
## IMN4 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## IN1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## IN2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## IN3 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## IPAN1.1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## IPAN1.2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## IPAN2.1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## IPAN2.2 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## IPAN3 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## IPAN4 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## MIX 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## Glia 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## SMC 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## sort_clusters
## preAnno2 13 16 23 6 19 24 3 20 15 21 25 22
## EMN1 0 0 0 0 0 33 0 0 0 0 0 0
## EMN2 0 0 0 0 0 0 0 0 0 0 0 0
## EMN3 0 0 0 0 0 0 0 0 0 0 0 0
## EMN4 0 0 0 0 0 0 0 0 0 0 0 0
## EMN5 0 0 0 0 0 0 0 0 0 0 0 0
## IMN1 9 0 0 0 0 0 0 0 0 0 0 0
## IMN2 404 22 0 0 0 0 0 0 0 0 0 0
## IMN3 0 306 42 0 0 0 0 0 0 0 0 0
## IMN4 0 0 212 0 0 0 0 0 0 0 1 0
## IN1 0 0 0 539 5 0 0 0 0 0 0 0
## IN2 0 0 0 7 301 0 0 0 0 0 0 0
## IN3 0 0 0 0 0 118 0 0 0 0 0 0
## IPAN1.1 0 0 0 0 0 0 611 1 0 0 0 0
## IPAN1.2 0 0 0 0 0 0 51 299 0 0 0 0
## IPAN2.1 0 0 0 0 0 0 0 1 294 5 0 0
## IPAN2.2 0 1 0 0 0 0 0 0 40 287 0 0
## IPAN3 0 0 1 0 0 0 0 0 0 0 95 0
## IPAN4 0 0 0 0 0 0 0 0 0 0 8 276
## MIX 0 0 0 0 0 0 0 0 0 0 0 0
## Glia 0 0 0 0 0 0 0 0 0 0 0 0
## SMC 0 0 0 0 0 0 0 0 0 0 0 0
markers.new.s <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
"Gfra2","Oprk1","Adamtsl1",
"Fbxw15","Fbxw24","Chrna7",
"Satb1","Cntnap5b",
"Gabrb1","Nxph1","Lama2","Lrrc7",
"Ryr3","Eda","Tac1",
"Kctd8","Ntrk2","Penk",
"Fut9","Nfatc1","Egfr","Mgll",
"Chrm3"
),
IMN=c("Nos1","Kcnab1",
"Gfra1","Etv1",
"Man1a","Airn",
"Adcy2","Cmah","Creb5","Vip","Pde1a",
"Ebf1","Gpc5"
),
IN=c("Npas3","Synpr","St18","Gal",
"Kcnk13",
"Neurod6","Moxd1","Sctr",
"Piezo1","Sst","Adamts9","Kcnn2",
"Calb2"),
IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
"Ngfr","Galr1","Il7","Aff2",
"Gpr149","Cdh6","Cdh8",
"Clstn2","Ano2","Ntrk3",
"Cpne4","Vwc2l","Cdh9","Scgn",
"Vcan","Cck","Piezo2","Kcnh7",
"Rerg","Bmpr1b","Skap1","Ntng1",
"Tafa2","Nxph2"),
CONTAM=c("Actl6b","Phox2b"),
pDEGs=c("Grid2","Ide","Btaf1","Cdh10","Slc15a2","Ccser1","Hdac9","Col25a1","Rspo2"))
pm.All_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "sort_clusters",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new.s
# Anno1
GEX.seur$Anno1 <- as.character(GEX.seur$seurat_clusters)
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(0,1,2,12)] <- "EMN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(11,17)] <- "EMN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(9)] <- "EMN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(10)] <- "EMN4"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(18)] <- "EMN5"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(4,5,14,8)] <- "IMN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(7,13)] <- "IMN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(16)] <- "IMN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(23)] <- "IMN4"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(6)] <- "IN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(19)] <- "IN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(24)] <- "IN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(3,20)] <- "IPAN1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(15,21)] <- "IPAN2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(25)] <- "IPAN3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(22)] <- "IPAN4"
GEX.seur$Anno1 <- factor(GEX.seur$Anno1,
levels = c(paste0("EMN",1:5),
paste0("IMN",1:4),
paste0("IN",1:3),
paste0("IPAN",1:4)
))
# Anno2
GEX.seur$Anno2 <- as.character(GEX.seur$seurat_clusters)
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(0,1,2,12)] <- "EMN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(11,17)] <- "EMN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(9)] <- "EMN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(10)] <- "EMN4"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(18)] <- "EMN5"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(4,5,14,8)] <- "IMN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(7,13)] <- "IMN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(16)] <- "IMN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(23)] <- "IMN4"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(6)] <- "IN1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(19)] <- "IN2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(24)] <- "IN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(3)] <- "IPAN1.1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(20)] <- "IPAN1.2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(15)] <- "IPAN2.1"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(21)] <- "IPAN2.2"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(25)] <- "IPAN3"
GEX.seur$Anno2[GEX.seur$Anno2 %in% c(22)] <- "IPAN4"
GEX.seur$Anno2 <- factor(GEX.seur$Anno2,
levels = c(paste0("EMN",1:5),
paste0("IMN",1:4),
paste0("IN",1:3),
paste0("IPAN",c(1.1,1.2,2.1,2.2,3,4))
))
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno1") +
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2")
markers.new.s <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
"Gfra2","Oprk1","Adamtsl1",
"Fbxw15","Fbxw24","Chrna7",
"Satb1","Cntnap5b",
"Gabrb1","Nxph1","Lama2","Lrrc7",
"Ryr3","Eda","Tac1",
"Kctd8","Ntrk2","Penk",
"Fut9","Nfatc1","Egfr","Mgll",
"Chrm3"
),
IMN=c("Nos1","Kcnab1",
"Gfra1","Etv1",
"Man1a","Airn",
"Adcy2","Cmah","Creb5","Vip","Pde1a",
"Ebf1","Gpc5"
),
IN=c("Npas3","Synpr","St18","Gal",
"Kcnk13",
"Neurod6","Moxd1","Sctr",
"Piezo1","Sst","Adamts9","Kcnn2",
"Calb2"),
IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
"Ngfr","Galr1","Il7","Aff2",
"Gpr149","Cdh6","Cdh8",
"Clstn2","Ano2","Ntrk3",
"Cpne4","Vwc2l","Cdh9","Scgn",
"Vcan","Cck","Piezo2","Kcnh7",
"Rerg","Bmpr1b","Skap1","Ntng1",
"Tafa2","Nxph2"),
CONTAM=c("Actl6b","Phox2b"),
pDEGs=c("Grid2","Ide","Btaf1","Cdh10","Slc15a2","Ccser1","Hdac9","Col25a1","Rspo2"))
pm.A1_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "Anno1",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All Anno1")
pm.A1_new.s
pm.A2_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "Anno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All Anno2")
pm.A2_new.s
I have noticed that annotations of EMN1-5 each time are a little
different from others ..
the boudaries are not that distinct to separate them,
better to do integration if have to measure the subtypes for multiple
datasets
GEX.seur$FB.info <- factor(as.character(GEX.seur$FB.info),
levels = c(paste0("CTL.",1:3),
paste0("CKO.",1:3)))
GEX.seur$rep <- paste0("rep",
gsub("CTL|CKO","",as.character(GEX.seur$FB.info)))
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno1") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno2") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.info,
clusters=GEX.seur$Anno1),
main = "Cell Count",
gaps_row = c(3),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.info,
clusters=GEX.seur$Anno1)),
main = "Cell Ratio",
gaps_row = c(3),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
stat_Anno1 <- GEX.seur@meta.data[,c("cnt","rep","Anno1")]
stat_Anno1.s <- stat_Anno1 %>%
group_by(cnt,rep,Anno1) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom.Anno1 <- stat_Anno1.s %>%
ggplot(aes(x = Anno1, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = c("skyblue","pink"), name="") +
labs(title="Cell Composition of SI data, Anno1", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=Anno1, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom.Anno1
glm.nb - anova.LRT
Sort.N <- c("CTL","CKO")
glm.nb_anova.Anno1 <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat_Anno1.s_N <- stat_Anno1.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat_Anno1.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat_Anno1.s_N$total <- total.N[paste0(stat_Anno1.s_N$cnt,
"_",
stat_Anno1.s_N$rep),"total"]
glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat_Anno1.s_N$Anno1)){
glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat_Anno1.s_N %>% filter(Anno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.Anno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}
glm.nb_anova.Anno1_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.Anno1)))
rownames(glm.nb_anova.Anno1_df) <- names(glm.nb_anova.Anno1)
colnames(glm.nb_anova.Anno1_df) <- gsub("X","C",colnames(glm.nb_anova.Anno1_df))
glm.nb_anova.Anno1_df
## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1
## CTLvsCKO 0.3141664 0.3103288 0.1787928 0.001014898 0.5353266 0.1853838
## IMN2 IMN3 IMN4 IN1 IN2 IN3
## CTLvsCKO 0.1070702 0.6577923 0.05522512 4.638542e-05 0.2279499 0.009895313
## IPAN1 IPAN2 IPAN3 IPAN4
## CTLvsCKO 0.4831668 0.8723714 0.2798343 0.8738828
round(glm.nb_anova.Anno1_df,5)
## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2 IMN3
## CTLvsCKO 0.31417 0.31033 0.17879 0.00101 0.53533 0.18538 0.10707 0.65779
## IMN4 IN1 IN2 IN3 IPAN1 IPAN2 IPAN3 IPAN4
## CTLvsCKO 0.05523 5e-05 0.22795 0.0099 0.48317 0.87237 0.27983 0.87388
cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.info,
clusters=GEX.seur$Anno2),
main = "Cell Count",
gaps_row = c(3),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.info,
clusters=GEX.seur$Anno2)),
main = "Cell Ratio",
gaps_row = c(3),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
stat_Anno2 <- GEX.seur@meta.data[,c("cnt","rep","Anno2")]
stat_Anno2.s <- stat_Anno2 %>%
group_by(cnt,rep,Anno2) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom.Anno2 <- stat_Anno2.s %>%
ggplot(aes(x = Anno2, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = c("skyblue","pink"), name="") +
labs(title="Cell Composition of SI data, Anno2", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=Anno2, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom.Anno2
glm.nb - anova.LRT
Sort.N <- c("CTL","CKO")
glm.nb_anova.Anno2 <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat_Anno2.s_N <- stat_Anno2.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat_Anno2.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat_Anno2.s_N$total <- total.N[paste0(stat_Anno2.s_N$cnt,
"_",
stat_Anno2.s_N$rep),"total"]
glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat_Anno2.s_N$Anno2)){
glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat_Anno2.s_N %>% filter(Anno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.Anno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}
glm.nb_anova.Anno2_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.Anno2)))
rownames(glm.nb_anova.Anno2_df) <- names(glm.nb_anova.Anno2)
colnames(glm.nb_anova.Anno2_df) <- gsub("X","C",colnames(glm.nb_anova.Anno2_df))
glm.nb_anova.Anno2_df
## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1
## CTLvsCKO 0.3141664 0.3103288 0.1787928 0.001014898 0.5353266 0.1853838
## IMN2 IMN3 IMN4 IN1 IN2 IN3
## CTLvsCKO 0.1070702 0.6577923 0.05522512 4.638542e-05 0.2279499 0.009895313
## IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2 IPAN3 IPAN4
## CTLvsCKO 0.5696465 0.3524691 0.8613843 0.6981254 0.2798343 0.8738828
round(glm.nb_anova.Anno2_df,5)
## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2 IMN3
## CTLvsCKO 0.31417 0.31033 0.17879 0.00101 0.53533 0.18538 0.10707 0.65779
## IMN4 IN1 IN2 IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2 IPAN3
## CTLvsCKO 0.05523 5e-05 0.22795 0.0099 0.56965 0.35247 0.86138 0.69813 0.27983
## IPAN4
## CTLvsCKO 0.87388
#saveRDS(GEX.seur,"./GEX0430.seur.pure_Anno.rds")